/*
 * Created on 2006-7-4
 */
package mathUtil;

/**
 * @author g0404403
 * @email wuw@nus.edu.sg The source code is copied from
 *        http://akiti.ca/Quad4Deg.html
 */
public class QuarticSolver {

    public QuarticSolver(double a, double b, double c, double d, double e) {
        this.solve(a, b, c, d, e);
    }

    public void solve(double a, double b, double c, double d, double e) {
        if (a == 0) {
            System.out
                    .println("The coefficient of the fourth power of x is 0. Please use the utility for a THIRD degree quadratic. No further action  taken.");
            return;
        } // End if a == 0

        if (e == 0) {
            System.out
                    .println("One root is 0. Now divide through by x and use the utility for a THIRD degree quadratic to solve the resulting   equation for the other three roots. No further action taken.");
            return;
        } // End if e == 0

        if (a != 1) {
            b /= a;
            c /= a;
            d /= a;
            e /= a;
        }

        double cb, cc, cd; // Coefficients for use with cubic solver
        double discrim, q, r, RRe, RIm, DRe, DIm, dum1, ERe, EIm, s, t, term1, r13, sqR, y1, z1Re, z1Im, z2Re;

        cb = -c;
        cc = -4.0 * e + d * b;
        cd = -(b * b * e + d * d) + 4.0 * c * e;

        if (cd == 0)
            System.out.println("cd = 0.");

        // Solve the resolvant cubic for y1
        q = (3.0 * cc - (cb * cb)) / 9.0;
        r = -(27.0 * cd) + cb * (9.0 * cc - 2.0 * (cb * cb));
        r /= 54.0;
        discrim = q * q * q + r * r;
        term1 = (cb / 3.0);

        if (discrim > 0) { // one root real, two are complex
            s = r + Math.sqrt(discrim);
            s = ((s < 0) ? -Math.pow(-s, (1.0 / 3.0)) : Math
                    .pow(s, (1.0 / 3.0)));
            t = r - Math.sqrt(discrim);
            t = ((t < 0) ? -Math.pow(-t, (1.0 / 3.0)) : Math
                    .pow(t, (1.0 / 3.0)));
            y1 = -term1 + s + t;
        } // End if (discrim > 0)
        else {
            if (discrim == 0) {
                r13 = ((r < 0) ? -Math.pow(-r, (1.0 / 3.0)) : Math.pow(r,
                        (1.0 / 3.0)));
                y1 = -term1 + 2.0 * r13;
            } // End if (discrim == 0)
            else { // else discrim < 0
                q = -q;
                dum1 = q * q * q;
                dum1 = Math.acos(r / Math.sqrt(dum1));
                r13 = 2.0 * Math.sqrt(q);
                y1 = -term1 + r13 * Math.cos(dum1 / 3.0);
            } // End discrim < 0
        } // End else discrim <= 0

        // At this point, we have determined y1, a real root of the resolvent
        // cubic.
        // Carry on to solve the original quartic equation

        term1 = b / 4.0;
        sqR = -c + term1 * b + y1; // R-squared

        RRe = RIm = DRe = DIm = ERe = EIm = z1Re = z1Im = z2Re = 0;

        if (sqR >= 0) {
            if (sqR == 0) {
                dum1 = -(4.0 * e) + y1 * y1;
                if (dum1 < 0) // D and E will be complex
                    z1Im = 2.0 * Math.sqrt(-dum1);
                else { // else (dum1 >= 0)
                    z1Re = 2.0 * Math.sqrt(dum1);
                    z2Re = -z1Re;
                }// End else (dum1 >= 0)
            } // End if (sqR == 0)
            else { // (sqR > 0)
                RRe = Math.sqrt(sqR);
                z1Re = -(8.0 * d + b * b * b) / 4.0 + b * c;
                z1Re /= RRe;
                z2Re = -z1Re;
            } // End else (sqR > 0)
        } // end if (sqR >= 0)
        else { // else (sqR < 0)
            RIm = Math.sqrt(-sqR);
            z1Im = -(8.0 * d + b * b * b) / 4.0 + b * c;
            z1Im /= RIm;
            z1Im = -z1Im;
        } // End else (sqR < 0)

        z1Re += -(2.0 * c + sqR) + 3.0 * b * term1;
        z2Re += -(2.0 * c + sqR) + 3.0 * b * term1;

        // At this point, z1 and z2 should be the terms under the square root
        // for D and E

        if (z1Im == 0) { // Both z1 and z2 real
            if (z1Re >= 0)
                DRe = Math.sqrt(z1Re);
            else
                DIm = Math.sqrt(-z1Re);
            if (z2Re >= 0)
                ERe = Math.sqrt(z2Re);
            else
                EIm = Math.sqrt(-z2Re);
        }// End if (zIm == 0)
        else { // else (zIm != 0); calculate root of a complex number********
            r = Math.sqrt(z1Re * z1Re + z1Im * z1Im); // Calculate r, the
                                                        // magnitude
            r = Math.sqrt(r);

            dum1 = Math.atan2(z1Im, z1Re); // Calculate the angle between the
                                            // two vectors
            dum1 /= 2; // Divide this angle by 2
            ERe = DRe = r * Math.cos(dum1); // Form the new complex value
            DIm = r * Math.sin(dum1);
            EIm = -DIm;
        } // End else (z1Im != 0)

        System.out.println();
        double x1Re = -term1 + (RRe + DRe) / 2;
        double x1Im = (RIm + DIm) / 2;
        double x2Re = -(term1 + DRe / 2) + RRe / 2;
        double x2Im = (-DIm + RIm) / 2;
        double x3Re = -(term1 + RRe / 2) + ERe / 2;
        double x3Im = (-RIm + EIm) / 2;
        double x4Re = -(term1 + (RRe + ERe) / 2);
        double x4Im = -(RIm + EIm) / 2;

        System.out.println("root1:" + x1Re + "," + x1Im);
        System.out.println("root2:" + x2Re + "," + x2Im);
        System.out.println("root3:" + x3Re + "," + x3Im);
        System.out.println("root4:" + x4Re + "," + x4Im);
    }

}
